%Version 1.0 Copyright William Devenport, Virginia Tech, 2/1/18
function CN ()
global ui flow;
close all
flow.g=1.4; %gamma
flow.xmax=0;
flow.pbpc=1;

ui.Figure=figure('Units','normalized','Name','Converging Nozzle','NumberTitle','off');
set(gcf,'Position',[0.25 0.25 0.40 0.5]);

ui.axes(1)=axes('Position',[0.1300    0.7093    0.6    0.2157],'Box','on');
ui.axes(2)=axes('Position',[0.1300    0.4096    0.6    0.2157],'Box','on');
ui.axes(3)=axes('Position',[0.1300    0.1100    0.6    0.2157],'Box','on');
ui.xlim=[-3 2];

ui.design=uicontrol(gcf,'Style','pushbutton','Units','normalized','Position',[0.77 0.12 0.2 0.05],'String','Design','Callback',@design);
ui.pbpcL=uicontrol(gcf,'Style','text','Units','normalized','Position',[0.775 0.90 0.06 0.05],'String','Pb/Pc','FontAngle','italic','FontSize',10);
ui.gammaL=uicontrol(gcf,'Style','text','Units','normalized','Position',[0.905 0.90 0.06 0.05],'String','Cp/Cv','FontAngle','italic','FontSize',10);
ui.pbpcV=uicontrol(gcf,'Style','text','Units','normalized','Position',[0.775 0.36 0.06 0.03],'String',num2str(flow.pbpc),'FontAngle','italic','FontSize',10);
ui.gammaV=uicontrol(gcf,'Style','text','Units','normalized','Position',[0.905 0.36 0.06 0.03],'String',num2str(flow.g),'FontAngle','italic','FontSize',10);
ui.pbpc=uicontrol(gcf,'Style','slider','Units','normalized','Position',[0.785 0.4 0.04 0.52],'Min',0.005,'Max',1,'SliderStep',[0.005 0.1]/0.995,'value',flow.pbpc,'Callback',@pbpcApply);
ui.gamma=uicontrol(gcf,'Style','slider','Units','normalized','Position',[0.915 0.4 0.04 0.52],'Min',1.2,'Max',1.7,'SliderStep',[0.01 0.1]/(1.7-1.2),'Value',flow.g,'Callback',@gammaApply);
ui.state=uicontrol(gcf,'Style','text','Units','normalized','Position',[0.18 0.93 0.5 0.05],'Fontsize',12);
flowplot();
end

function flowplot ()
global ui flow;

g=flow.g;
xmax=flow.xmax;
pbpc=flow.pbpc;

%Nozzle shape
x=-3:.1:0;
y=1+x.^2/4.5;

%Design exit Mach and pressure
meD=1;
pbpcD=pp0(meD,g);flow.pbpcD=pbpcD;
%Subsonic choked exit Mach and pressure
meC=1;
pbpcC=pbpcD;flow.pbpcC=pbpcC;

axes(ui.axes(1)); % plot nozzle
ymax=max([3 y(end)*1.5]);
fill([x x(end) x(1)],[y  ymax ymax],[0.85 0.85 0.85]);hold on
plot(x,y,'k');ylim([0 ymax]);ylabel('A/A_t');xlim(ui.xlim);
plot([x(end) x(end)],[y(end) ymax],'k');hold off;

if pbpc>pbpcC %subsonic
    mexit=m_pp0(pbpc,g);
    aeas=aas(mexit,g);
    m=m_aas(y*aeas,g,0);m(find(y*aeas<=1))=1;
    ppc=pp0(m,g);
    axes(ui.axes(2));plot(x,m,'b');ylabel('M');yl=ylim;ylim([0 yl(2)]);xlim(ui.xlim);
    axes(ui.axes(3));plot(x,ppc,'b');ylabel('p/p_c');yl=ylim;ylim([0 yl(2)]);xlim(ui.xlim);
    hold on;plot([x(end) ui.xlim(2)-0.5],[ppc(end) pbpc],'b:',[ui.xlim(2)-0.5 ui.xlim(2)],[pbpc pbpc],'b');hold off;
    if pbpc<1
        axes(ui.axes(1));hold on; %plot flow
        hflow=fill([x ui.xlim(2) ui.xlim(2) ui.xlim(1)],[y y(end) 0 0],[0.85 0.85 1]);set(hflow,'Linestyle','none');
        hold off 
        set(ui.state,'String','Subsonic Flow');
    else
        set(ui.state,'String','No Flow');
    end
elseif pbpc<=pbpcC %underexpanded and design
    up=find(x<0);dn=find(x>=0);
    m=[m_aas(y(up),g,0) m_aas(y(dn),g,1)];
    ppc=pp0(m,g);
    axes(ui.axes(2));plot(x,m,'b');ylabel('M');yl=ylim;ylim([0 yl(2)]);xlim(ui.xlim);
    axes(ui.axes(3));plot(x,ppc,'b');ylabel('p/p_c');yl=ylim;ylim([0 yl(2)]);xlim(ui.xlim);
    hold on;plot([x(end) ui.xlim(2)-0.5],[ppc(end) pbpc],'b:',[ui.xlim(2)-0.5 ui.xlim(2)],[pbpc pbpc],'b');hold off;
    if pbpc==pbpcD
        set(ui.state,'String','Design Condition, subsonic choked');
        axes(ui.axes(1));hold on; %plot flow
        hflow=fill([x ui.xlim(2) ui.xlim(2) ui.xlim(1)],[y y(end) 0 0],[0.85 0.85 1]);set(hflow,'Linestyle','none');
        hold off 
    elseif pbpc<pbpcD
        set(ui.state,'String','Underexpanded Flow');
        m1=meD;m2=m_pp0(pbpc,g);theta=nu(m2,g)-nu(m1,g);mu1=asin(1/m1);mu2=asin(1/m2)-theta;
        axes(ui.axes(1));hold on; %plot flow and waves
        hflow=fill([x ui.xlim(2) ui.xlim(2) ui.xlim(1)],[y y(end)+(ui.xlim(2)-x(end))*tan(theta) 0 0],[0.85 0.85 1]);set(hflow,'Linestyle','none');
        xs=[x(end) x(end)+y(end)*cot(mu1) x(end)+y(end)*cot(mu2)];ys=[y(end) 0 0];
        if mu2<0
            xs=[x(end) ui.xlim(2) ui.xlim(2)];ys=[y(end) y(end)-(ui.xlim(2)-x(end))*tan(mu1) y(end)-(ui.xlim(2)-x(end))*tan(mu2)];    
        end
        hwave=fill(xs,ys,[0.85-0.85*tanh(m2-m1) 0.85 1-tanh(m2-m1)]);set(hwave,'Linestyle','none');
       hold off;   
    end    
end
end

function pbpcApply (object, eventdata)
global ui flow;
flow.pbpc=get(ui.pbpc,'Value');
flowplot();
set(ui.pbpcV,'String',num2str(flow.pbpc));
set(ui.pbpc,'Value',flow.pbpc);
end

function gammaApply (object, eventdata)
global ui flow;
flow.g=get(ui.gamma,'Value');
flowplot();
set(ui.gammaV,'String',num2str(flow.g));
set(ui.gamma,'Value',flow.g);
end

function design (object, eventdata)
global ui flow;
flow.pbpc=flow.pbpcD;
flowplot();
set(ui.pbpcV,'String',num2str(flow.pbpc));
set(ui.pbpc,'Value',flow.pbpc);
end

function r=aas(m,g)
r=1./m.*(2/(g+1)*(1+(g-1)/2*m.^2)).^((g+1)/2/(g-1));
end

function r=pp0(m,g)
r=(1+(g-1)/2*m.^2).^(-g/(g-1));
end

function r=p02p01(m,g)
r=((g+1)*m.^2./((g-1)*m.^2+2))^(g/(g-1)).*((g+1)./(2*g*m.^2-(g-1))).^(1/(g-1));
end

function r=p2p1(m,g)
r=1+2*g/(g+1)*(m.^2-1);
end

function n=nu(m,g)
   n=sqrt((g+1)/(g-1))*atan(sqrt((g-1)/(g+1)*(m.^2-1)))-atan(sqrt(m.^2-1));
end

function r=m_pp0(pp0,g)
r=sqrt((pp0.^(-(g-1)/g)-1)*2/(g-1));
end

function m=m_aas(a,g,super)  %super=1 returns supersonic solutions, else subsonic
m=zeros(size(a));
m(find(a<1))=nan;
m(find(a==1))=1;

if super==1
    minit=1e5;
else
    minit=1e-5;
end

for n=1:length(a(:))
    if m(n)==0
        m0=1;m1=minit;
        while(abs(m1-m0)>1e-6)
            m0=m1;
            f=aas(m0,g)-a(n);
            fp=(2/(g+1)*(1+(g-1)/2*m0^2)).^((g+1)/2/(g-1)-1)-aas(m0,g)/m0;
            m1=m0-f/fp;
            if m1<0
                m1=m0/2;
            end
        end
        m(n)=m1;
    end
end
end




